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USE  OF  THE  TENSOR  PRODUCT  FOR  NUMERICAL  WEATHER 
PREDICTION  BY  THE  FINITE  ELEMENT  METHOD  -  PART  2. 


Introduction 

This  is  the  second  installment  of  a  report-pair  concern¬ 
ing  implementation  of  tensor  product  factoring  of  coeffi¬ 
cient  matrices  in  applications  of  the  finite  element  method 
to  numerical  weather  prediction.  It  was  noted  in  Part  1 
(Ref.  1)  that  these  techniques  were  introduced  in  numerical 
weather  prediction  by  Staniforth  and  Mitchell  (Ref.  2). 
Discussed  in  Part  1  are  applications  in  which  the  "mass" 
matrix  for  a  grid  such  as  that  shown  in  Fig.  1  is  factored 
as  the  tensor  product  of  two  matrices. 


n  columns 

Fig.  1.  Node  numbering  and  spacing. 

One  of  these  matrices  (MA)  depends  solely  upon  the  nodal- 
spacing  in  the  east-west  direction  (a^ )  and  the  other  (MB) 
depends  only  on  the  north-south  spacing  (b^).  We  began  with 
the  set  of  simultaneous  linear  equations 

M  w  =  v,  <1> 

where  M  (the  "mass"  matrix)  is  symmetric,  ne  x  ne,  and  w  and 
v  are  column  vectors  of  height  ne.  M  and  v  are  input  quan¬ 
tities  and  w  is  sought.  The  tensor  product  representation 
of  M  is 

M  =  MB  *  MA,  <2> 

where  MB  and  MA  are  tridiagonal,  symmetric  matrices,  e  x  e 
and  n  x  n,  respectively.  (The  tensor  product  and  matrices 
MA  and  MB  are  defined  in  Appendix  A. )  This  representation 
allowed  <1>  to  be  rewritten  as 

MA  W  MB  =  V,  <3> 

where  W  is  n  x  e  and  the  successive  columns  are  subvectors 


of  w  corresponding  to  the  rows  of  Fig.  1.  V  is  also  n  x  e 
and  similarly  derived  from  v.  Boundary  conditions  consid¬ 
ered  were  a  cyclic  condition  in  the  east-west  direction  and 
either  homogeneous  Neumann  conditions  (normal  derivative 
zero)  or  nonhomogeneous  Dirichlet  conditions  (specified 
nonzero  values)  on  the  northern  and  southern  edges. 

The  present  report  discards  the  cyclic  east-west  boundary 
condition  and  deals  with  two  cases: 

(1)  Solutions  of  <3>  with  nonhomogeneous  Dirichlet  con¬ 
ditions  on  all  four  edges; 

(2)  Solution  of  Poisson's  equation  for  the  same  region 
with  nonhomogeneous  Dirichlet  conditions  on  all  four 
edges . 

Mass  Matrix  ^  Dirichlet  Boundary  Conditions 

Effects  of  the  Dirichlet  boundary  conditions  on  the  solu¬ 
tion  process  are  most  readily  understood  by  considering  the 
following  partitioned  form  of  <1>: 


In  <4>  the  w  vector  has  been  rearranged  so  that  all  of  the 
boundary  values  are  in  the  subvector  W|j  and  the  interior 
("center")  values  are  in  wc .  A  similar  reordering  has  been 
applied  to  v  and  M.  If  the  boundary  values  of  w  are  pre¬ 
scribed,  then  Wjj  is  known  and  only  wc  remains  to  be  found. 
Expanding  the  lower  partition  of  <4>  and  placing  the  known 
terms  on  the  right  gives 

M22WC  =  VC  ■  «2lWb»  <5> 
or,  letting  v  '  =  vc  *  Mziw^,  we  have 

M22wc  =  vc' .  <5’> 

We  consider  now  how  the  strategy  just  described  can  be 
applied  when  the  tensor  product  resolution  of  M  has  been 
used  to  convert  <1>  into  <3>.  In  the  matrix  W  the  pre¬ 
scribed  boundary  values  occupy  the  first  and  last  columns 
and  the  top  and  bottom  rows.  Denote  this  border  matrix, 
including  an  (n-2)  x  (e-2)  null  matrix  inside,  by  WB. 


Calculate 


VB  =  MA  WB  MB 


<6> 


and  now  form 


V'  =  V  -  VB. 


<7> 


Now  define  a  set  of  submatrices  MAI,  MB1,  Wl,  and  VI 
obtained  from  MA,  MB,  W,  and  V' ,  respectively,  by  removing 
the  first  and  last  columns  and  the  top  and  bottom  rows.  The 
reduced  problem  becomes 

MAI  Wl  MB1  =  VI  <8> 


As  described  in  Ref.  1,  <8>  may  be  solved  by  standard  Gaus¬ 
sian  elimination  procedures.  A  computer  program  (GAUSS4) 
which  carries  out  these  calculations  is  listed  in  Appendix 
B.  The  subroutines  of  GAUSS4  are  designed  for  substitution 
in  the  program  devised  by  Hinsman  (Ref.  5). 


Poisson’ s  Equation  ^  Dirichlet  Boundary  Conditions 

As  noted  above,  Staniforth  and  Mitchell  (Ref.  2)  appear 
to  have  been  first  in  applying  the  tensor  product  resolution 
to  Poisson's  equation  in  a  numerical  weather  prediction 
problem  using  the  finite  element  method.  Additional  detail 
is  given  in  earlier  papers  by  Dorr  (Ref.  3)  and  by  Lynch, 
Rice,  and  Thomas  (Ref.  4). 

Finite  element  discretization  of  Poisson's  equation  for 
the  region  of  Fig.  1  results  in  a  set  of  simultaneous  linear 
equations  which  may  be  written  in  matrix  form  as 

K  w  =  v,  <9> 

where  vectors  v  and  w  are,  respectively,  given  and  unknown. 
As  for  <1>,  each  has  length  ne  and  the  coefficient  matrix  K 
is  ne  x  ne ,  symmetric,  sparse,  and  block  tridiagonal.  K  is 
called  the  "stiffness"  matrix  in  finite  element  parlance. 

It  is  easily  shown  that  K  is  expressible  as  the  sum  of 
two  tensor  products  as  follows: 

K  =  MB  *  SA  +  SB  *  MA.  <10> 

The  new  matrices  SA  and  SB  are  symmetric,  tridiagonal  and 
depend  only  on  the  a^  and  b ,  respectively.  Explicit  formu¬ 
las  for  SA  and  SB  are  given  in  Appendix  A. 


Using  the  definition  of  the  tensor  product  and  again  con¬ 
verting  the  vectors  w  and  v  into  the  n  x  e  rectangular 
matrices  W  and  V,  <9>  may  be  written  as 

SA  W  MB  +  MA  W  SB  =  V.  <11> 

Before  solving  <11>  we  must  first  take  account  of  the  Diri- 
chlet  boundary  conditions  on  the  four  edges  of  the  region. 
As  in  solving  <3>,  the  given  boundary  values  are  in  the 
first  and  last  columns  and  top  and  bottom  rows  of  W.  As 
before,  we  let  WB  be  an  n  x  e  matrix  containing  the  given 
boundary  values,  together  with  zeros  at  locations  corre¬ 
sponding  to  interior  nodes.  Calculate 

VB  =  SA  WB  MB  ♦  MA  WB  SB,  <12> 

and  then  form 


V'  =  V  -  VB.  <7> 

The  remaining  step  again  parallels  that  used  when  applying 
the  Dirichlet  boundary  conditions  to  <3>.  Specifically,  we 
introduce  submatrices  MAI,  MB1 ,  SAl,  SBl,  Wl,  and  VI 
obtained  from  MA,  MB,  SA,  SB,  W,  and  V’,  respectively,  by 
removing  the  first  and  last  columns  and  the  top  and  bottom 
rows.  The  reduced  problem  becomes 

SAl  Wl  MB1  +  MAI  Wl  SBl  =  VI.  <13> 


To  solve  <13>  we  first  need  the  complete  solution  of  the 
eigenproblem 

SBl  p.j  =  A.j  MB1  p^  ,  <14> 

where  p^  is  the  ith  eigenvector  and  A^  is  the  corresponding 
eigenvalue.  We  write  the  complete  solution  in  the  form 

SBl  P  =  MB1  P  A,  <14’> 

where  P  is  the  (e-2)  x  (e-2)  modal  matrix  whose  columns  are 
the  p..  and  A  is  the  (diagonal)  spectral  matrix  whose  ele¬ 
ments  are  the  Aj  .  We  specify  that  the  modal  matrix  is  nor¬ 
malized  so  that 

PT  MBl  P  =  I,  <15> 

where  I  is  the  identity  matrix  of  order  e-2  and  PT  is  the 
transpose  of  P.  If  both  sides  of  <13>  are  postmultiplied  by 
P  and  <14’>  is  used  to  replace  SBl  P,  <13>  becomes 

SAl  Wl  MBl  P  ♦  MAI  Wl  MBl  PA  =  VI  P.  <16> 

Let  X  =  Wl  MBl  P  and  U  =  VI  P,  then  <16>  is  equivalent  to 
(SAl  ♦  A.j  MAI)  x^  =  u..  ,  i  =  1,  e-2,  <17> 
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where  and  are,  respectively,  the  ith  columns  of  X  and 
U.  Since  the  coefficient  matrix  in  <17>  is  tridiagonal,  the 
Gaussian  elimination  process,  i.e.,  factoring,  forward 
reduction,  and  back- substitution,  is  computationally  econom¬ 
ical.  The  final  step  consists  of  a  matrix  multiplication  to 
obtain 

W1  =  X  PT.  <18> 

Since  W1  contains  the  w  values  at  all  interior  nodes  and  the 
boundary  values  were  known  in  advance,  the  solution  is  com¬ 
plete.  A  FORTRAN  program  (GAUSS5)  which  implements  the  ten¬ 
sor  product  solution  for  Poisson's  equation  is  given  in 
Appendix  C. 

Operation  Counts  and  Storage  Requirements  ^  Poisson's  Equa¬ 
tion 

In  Ref.  1  comparisons  of  floating  point  operation  counts 
and  storage  requirements  were  made  for  solutions  of  <1>. 
Substitution  of  the  boundary  conditions  considered  here  in 
place  of  those  considered  in  Ref.  1  has  a  negligible  effect 
on  both  operation  counts  and  storage  requirements.  Accord¬ 
ingly,  no  further  comparison  is  given  here  for  solutions  of 
<1>. 

Solution  of  Poisson's  equation  (<9>)  using  the  tensor 
product  resolution  <10>  of  K  is  more  costly  in  terms  of  com¬ 
putation  and  storage  than  the  previously  studied  applica¬ 
tions  to  <1>.  In  Table  1  the  number  of  floating  point  oper¬ 
ations  and  the  required  number  of  coefficient  matrix  storage 
locations  are  compared  for  three  different  solution  methods. 
These  are  SOR  (successive  over- relaxation) ,  SKY  (skyline 
storage  and  Gauss  elimination),  and  TENSOR  (the  scheme 
described  above).  A  floating  point  operation  is  defined  to 
be  one  multiplication  (or  division)  plus  one  addition  (or 
subtraction).  The  exact  operation  counts  would  be  polynomi¬ 
als  in  n  and  e.  Only  the  highest  degree  terms  are  given  in 
the  table.  Since  it  is  not  possible  to  predict  the  number 
of  iterations  per  solution  using  SOR,  the  operation  count 
given  for  that  algorithm  is  for  a  single  iteration.  In 
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Table  1  a  storage  location  corresponds  to  8  bytes.  For  the 
comparison  it  is  assumed  that  each  floating  point  number 
requires  8  bytes  of  storage  and  an  integer  requires  4  bytes. 
The  storage  requirement  given  for  SOR  is  based  on  the  com¬ 
pact  storage  scheme  described  by  Franke  and  Salinas 
(Ref.  6). 


TABLE  1.  Operation  Counts  and  Storage  Requirements. 


Note:  1.  Number  of  operations  per  iteration. 


It  is  perhaps  surprising  to  note  that  the  number  of  oper¬ 
ations  for  TENSOR  is  no  fewer  than  for  SKY.  Turning  atten¬ 
tion  to  storage  requirements  reveals  that  for  a  large  prob¬ 
lem  (e  =  n  =  100,  say)  the  SKY  storage  requirement  for  the 
stiffness  matrix  is  8  megabytes,  compared  with  1  megabyte 
for  SOR  and  80  kilobytes  for  TENSOR.  It  is  this  comparison 
which  is  the  compelling  reason  for  preferring  TENSOR.  It  is 
acknowledged  that  there  is  overhead  associated  with  the  one¬ 
time  solution  of  the  eigenvalue  problem  <14>,  but  the  tri¬ 
diagonal  form  of  matrices  SB1  and  MB1  makes  the  amount  of 
computation  comparable  with  that  required  for  a  single  solu¬ 
tion  of  Poisson's  equation.  Since  two  solutions  of  Pois¬ 
son's  equation  are  required  at  each  time  step,  the  overhead 
is  clearly  negligible. 

It  is  not  feasible  to  make  a  definitive  comparison 
between  the  number  of  operations  required  for  SOR  and  those 
required  for  the  other  two  algorithms.  If  the  number  of 
iterations  is  less  than  0.2  e,  then  SOR  will  be  more  econom¬ 
ical  and  the  storage  tradeoff  would  need  to  be  weighed. 


Conclusions 


It  has  been  demonstrated  that  Dirichlet  boundary  condi¬ 
tions  on  all  edges  of  the  region  are  easily  incorporated  in 
solution  processes  which  use  tensor  product  resolution  of 
the  coefficient  matrix.  For  very  large  problems  the  tensor 
product  algorithm  uses  much  less  core  storage  than  alterna¬ 
tive  choices.  The  computational  expense  of  a  solution  to 
Poisson's  equation  is  substantially  the  same  for  Gaussian 
elimination  and  for  the  tensor  product  scheme.  It  is 
expected  that  successive  over- relaxation  is  almost  always 
more  expensive. 
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APPENDIX  A  -  TENSOR  PRODUCT  AND  MATRIX  DEFINITIONS 

Tensor  Product  The  tensor  product  of  matrices  C  and  D 

may  be  represented  in  block  partition  form  as 


"ciiD 

Cl2  0 

Cl 3  0‘ 

C  *  D 

=  c2i  D 

c22  D 

C  2  3  0 

,c3i  D 

c32  0 

c33  D, 

where  the  c^j  are 

the  elements  of 

C. 

Note  that, 

if  C  and  D 

have  dimensions  r 

x  s  and 

rt 

X 

c 

respectively , 

the  tensor 

product  has  dimensions  rt  x  su. 

Definitions  for  matrices  MA  and  SA  are  given  below.  The 
corresponding  expressions  for  MB  and  SB  may  be  obtained  by 
substituting  "b"  for  "a"  throughout  and  replacing  n  by  e. 
(Symbols  aj  and  bj  are  defined  in  Fig.  1.) 

2a i  3i  0  0 

1  ai  2(ai+a2)  a2  0 

m  =  6  0  a2  2(a2+a3)  a3 

(n=4)  _  0  0  a3  2a3  _ 


SA  = 
(n=4) 


I  +  I-I 

di  32  a2 

.1  i+I 

312  »3 
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APPENDIX  B 


PROGRAM  to  SOLVE  M  w  =  v  with  DIRICHLET  BOUNDARY  CONDITIONS 
Listing:  GAUSS4  FORTRAN 

C  MAIN  PROGRAM  MASS  MATRIX  USING  TENSOR  PRODUCT  RESOLUTION 
C 

C  THIS  PROGRAM  IS  DESIGNED  TO  TEST  THE  SCHEME  (TENSOR) 

C  WHICH  RESOLVES  THE  MASS  MATRIX  INTO  A  TENSOR  PRODUCT  IN 
C  ORDER  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS  M  w  =  v.IN 
C  THIS  PROGRAM  THERE  ARE  DIRICHLET  BOUNDARY  CONDITIONS  ON 
C  ALL  4  EDGES  OF  THE  REGION.  THE  PRESCRIBED  BOUNDARY 
C  VALUES  ARE  GIVEN  IN  THE  CORRESPONDING  LOCATIONS  IN  V. 

C  THE  SUBROUTINES  MAY  BE  INSERTED  IN  THE  PROGRAM  DEVISED 
C  BY  HINSMAN . 


IMPLICIT  REAL*8 (A-H ,0- Z ) 
COMMON/ CM 1A/ NLAT, NLONG 


COMMON/ CM8 
COMMON  AG? 
DIMENSION 
READ ( 5,  *) 

latx=n£at* 
WRITE (6 .10 


NLONG,  NLAT 


iENSION  V(ZP) 

D(5,  * ) NLONG ,  NLAT 
X=NLAT* I 
TEX6.10Q0) 

MAT ?  /  /  /  MASS  MATR 
TE£  6  ^ 100 1 ) NLONG , NLAT 


AD(ZK) ,GBD(ZL) ,MA(ZM) ,MB(ZN) 


MASS  MATRIX  -  TENSOR  PRODUCT  RESOLUTION’ 


1000 


503  FORMAT?/,;  B:  , (24F3.0)) 

500  FORMAT ?  /  j  A:  ’  (24F3.0J) 

1001  FORMAT ? ’  NLONG  =, 13,'  NLAT=',I3,/) 

C  CONSTRUCT  FACTORS,  GAD  AND  GBD,  OF  MASS  MATRIX 
CALL  AMTRX3 


1002 

1004 


CALL  AMTRX3 
WRITEf6.501)AG 
FORMAT?/,’  '  AG 
WRITE? 6 , $04 )BG 
FORMAT?//  BG 
WRITE? 6 , 1002  )GA 
FORMAT?//  GAD 
WRITE?6 , 1004 )GB 
FORMAT?///  GB 
WRITE ? 6, 1003 )MA 
WRITE? 6 , 10Q6JMB 

v - 7mt  at:  1  1  *MT  n\T 


(///AG:  ,  ( 12F4. 1) ) 

6.$04)BG 

(//  BG:  ,  ( 12F4. 1)  ) 

te&/'c3x,6F7-3>) 

^C)SD''/-(3x-6F7-3)) 


K= (NLAT- 1)*NL 
L=NLAT*NLONG 
READ £5  * TV 
WRITE? 6, 5 10 )V 
CORRECT  RIGHT -HA 
LONGM= NLONG- 1 


*  )V 

0,5lO)V 

ight-ha: 


ND  SIDE  FOR  DIRICHLET  CONDITION 


J>?2* 
CU=GBD( 
CL=GBD( 
GBD? 3?= 
>p*L 


-  1)*V( J- l)+GAD(2*J-2) 

J - 1 ) +  GAD ( 2* J - 2 ) *V ( L+ J ) 


S^LATX-l) 
LATX- 1)  =  0 . 
)£ftLONG+2)= 


m 


2*GAD?2*NL0NG-1) 

GBD?  3  j  =  CU 
GBD£2*LATX-1)=CL 
WRITE (6, 5 10 )V 
PERFORM  LDLT  FACTORING 
CALL  FACT1 ( GAD , NLONC 
CALL  FACT 1? GBD ! LATX, 
WRITE (6, 1002 TgAd 
WRITE (6  1004) GBD 


ins(j-2> 


CALL  FACT 1( GAD, 
CALL  FACT 1 (GBD ! 
WRITE (6, 1002 TgA 
WRITE  (6  J  -  '*'*** 


)F  GAD  AND  GBD 
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( 


1003 

1006 


PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
FACTORS  OF  GAD 

CALL  BACKA1(GAD,V) 

WRITE(6,510}v 

PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
FACTORS  OF  GBD 

CALL  B ACKB 1 ( GBD , V ) 

WRITEC6 ,510jV 

FORMAT!/,;  V:  ’ , 5F8 . 2./ . (4X , 5F8 . 2 j j 
3  FORMAT!/  MA:Y,2X,m3) 

6  FORMAT!/  ’  MB:  ,2X,36I3) 

STOP 


SUBROUTINE  FACT1(A,NN) 

C  'SUBROUTINE  FACTi  PERFORMS  L*D*LT ’ FACTORING ’ ON  A  SUBMATRIX 
C  OF  A  SYMMETRIC  TRIDIAGONAL  MATRIX  STORED  IN  SKYLINE  FORM. 
C  THE  SUBMATRIX  IS  FORMED  BY  OMITTING  THE  FIRST  AND  LAST 
C  COLUMNS  AND  ROWS  OF  THE  INPUT  MATRIX. 

C  .-  -  INPUT  VARIABLES  -  - 

C  .  A(NWK)  =  INPUT  MATRIX  STORED  IN  COMPACTED  FORM 

C  .  NN  =  NUMBER  OF  COLUMNS  (OR  ROWSJ  IN  INPUT  MATRIX  . 

C  .  NWK  =  NUMBER  OF  ELEMENTS  BELOW  SKYLINE  (2*NN  -  1)  . 

C  .  -  -  OUTPUT  -  -  . 

C  .  A (NWK)  =  D  AND  L  -  FACTORS  OF  INPUT  SUBMATRIX 


AND  L  -  FACTORS  OF  INPUT  SUBMATRIX 


IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  A?l) 


PERFORM  L*D*LT  FACTORIZATION  OF  STIFFNESS  MATRIX 
LONGMM=NN-2 

TEM?=Af  2^+l?/A?2*(J  -  1 ) ) 

WRITE ( IOUT , 2000 )N | A (KN ) 

TEMP 

*  c  m 


50 

2000 


-A  (2^. 
A(KN) 


2E20.12) 

RETURN 


/.T  STOP  -  MATRIX  NOT  POSITIVE  DEFINITE'  // 
ITIVE  PIVOT  FOR  EQUATION  ,14,///  PIVOT  = 


*S*S*^iX*-iXi*T(f************-**'********  •*■***■*■*** 

C  THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 
C  SUBSTITUTION  USING  THE  FACTORS  OF  GAD 


IMPLICIT 
COMMON/C 
DIMENS 10 


IT  REAL*8 (A-H.O-Z ) 


DEFINE  LIMITS  FOR  DO -LOOPS 

NTM=NLAT- 1 
LONGM=NLONG- 1 
LONGMM=  NLONG - 2 

REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 

DO  100  K= 1 , NTM 
DO  20  J=3,f.ONGM 

V(K%LONG+J)=V(K*NLONG+J)-V(K*NLONG+J-l)*A(2*J-l) 
DIVIDE  BY  DIAGONAL  ELEMENTS 


DO  40  J= 1 , LONG MM 
V(K*NLONG+ J* 1 )=V 


V (K*NLONG+ J+ 1) /A(2*J ) 


BACK- SUBSTITUTE 
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DO  60  J=3 .LONGM 
L=f^+lJ*NtONG-J+l 
M=  2*  ( NLONG -  J  j  +  3 

60  v(l)=v7l)-v(l+i)*a(m) 

100  CONTINUE 

RETURN 
END 
C 

Q  •hirirlflrkititieifk-b'iritififkiciS'Jfkii-Ifk'ititifie-irSfltifk'ieieirirkickififirleifitiirle'U-ie-le'kt'c'ict'c'Jcie-f: 

c 
c 
c 
c 


c 

c 

c 


c 

c 

c 


*'******3r**"w5**'****Vr*'Sw*i5r  #*#*******************■***•*'****#'*** 

THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 
SUBSTITUTION  USING  THE  FACTORS  OF  GBD. 

IMPLICIT  REAL*8 (A-H ,0- Z ) 

COMMON/ CM 1A/ NLAT . NLONG 
DIMENSION  A(1),V(1) 

DEFINE  NEEDED  INDEX  VARIABLES 

LATX=NLAT+ 1 
LONGM= NLONG -1 

REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 


20 

C 

C 


DO  100^K=  2^L0NGM 

Jills’- 1  j  ^NLONG )  =  V  ( K+  (  J  - 1 )  *NLONG  ) - V (K+ ( J - 2 )  *NLONG ) 


DIVIDE  BY  DIAGONAL  ELEMENTS 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 


2 

C 


THIS  SUBROUTINE  FORMS  THE  MASS  MATRIX  IN  THE  FORM  OF  A 
TENSOR  PRODUCT  OF  THE  GBD  MATRIX  AND  THE  GAD  MATRIX. 

THE  FIRST  OF  THESE  IS  NLAT  +  1  BY  NLAT  ♦  1.  SYMMETRIC. 

AND  TRIDIAGONAL.  THE  SECOND  IS  NLONG  BY  NLONG.  SYMMET¬ 
RIC.  AND  TRIDIAGONAL.  NOTE  THAT  THERE  ISNO  CYCLIC 
BOUNDARY  CONDITION  IN  THE  EAST- WEST  DIRECTION.  BOTH  GBD 
AND  GAD  ARE  STORED  IN  SKYLINE  VECTOR  FORM  (UPPER  TRIANGLE 
WITH  SPACE  FOR  FILL-IN).  INTEGER  ADDRESS  VECTORS  MB  AND 
MA  ARE  ALSO  GENERATED . 


MB(ZN) 


LATX=NLAT+ 1 
LONGM= NLONG -1 

FIND  BG  =  (ELEMENT  HEIGHT)/ 6 
LONGM= NLONG -1 
DO  2  J=1.NLAT 

Bgi  ^!BIkL0NGM*  (  J  -  1 )  )  /  3  * 

GBD(l7=2.*BGQ) 

DO  4  J=2.NLAT 
K=  2* ( J - 1 ) 

GBD(K)=2:*(BG(J-1)+BG(J)) 
GBDlKtl?  “***  " 


GBD (2* 
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10 

c 


12 


16 


18 


GBD 

FIND 

DO 


GAD 

GAD 

GAD 

GAD 


K)=2^*(AG(J-1)+AG(J) ) 

1?MGM)  =  ?*AG(L0NGM) 
IRY  V 


_ ^2*L0NGM+ 1 

GENERATE  direct 
MB ( 1 1  =  1 
DO  16 
MB 


longm; 

’ECTOR^ 


__  =1,NLAT 

T+  1  ^  ^  0  n'  T 

N1AT+ 2 ) =2* (NLAT+ 1 ) 

18  J=2,NL0NG 
j\s2*(J- 1) 

,ONG+l)  =  2*NLONG 


MA(j) 
MAINL- 
RETURN 
END 
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APPENDIX  C 

PROGRAM  -  POISSON'S  EQUATION  with  DIRICHLET  BOUNDARY 

CONDITIONS 


Listing 


GAUSS 5  FORTRAN 


MAIN  PROGRAM  STIFFNESS  MATRIX  USING  TENSOR  PRODUCT 
RESOLUTION 

THIS  PROGRAM  IS  DESIGNED  TO  TEST  THE  SCHEME  WHICH 
RESOLVES  THE  STIFFNESS  MATRIX  INTO  A  SUM  OF  TWO  TENSOR 
PRODUCTS  IN  ORDER  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS 
K  W  =  V.  THERE  ARE  DIRICHLET  BOUNDARYCONDITIONS  ON 
ALL  4  EDGES  OF  THE  REGION.  THE  PRESCRIBED  BOUNDARY 
VALUES  ARE  GIVEN  IN  THE  CORRESPONDING  LOCATIONS  IN  V. 
THE  SUBROUTINES  MAY  BE  INSERTED  IN  THE  PROGRAM  DEVISED 
BY  HINSMAN. 


IMPLICIT  REAL*8 (A-H , 0- Z ) 


V/NLA 

'AfZl 


WVI  «  1VH  /  W*  l  a;  |V  l 


SBI(ZL) 


503  FORMA' 
500  FORMA' 


!mo9°) 

§  503)B. 


STIFFNESS  MATRIX  -  TENSOR  PRODUCT 
LONG , NLAT 


LATX=NLAT+1 
WRITE? ©, 1000 ) 

1000  FORMAT (//.  ;  STIFFNESS  MATRIX  -  TENSOR  PRODUCT 

lRESOLUTIO*'  / ) 

WRITE? 6 , 100 1) NLONG , NLAT 

read?5,Ma.b 

WRITE? 6, SOdj A 
WRITE? 6’ 503 )B 

503  FORMAT?/,;  B:  , (24F3.0)) 

500  FORMAT?/,'  A:  ',?24F3,0)) 

1001  FORMAT?’  NLONG  =  \  13  /  NLAT  =’  13./) 

C  CONSTRUCT  FACTORS,  GAl,  GB1,  SA1 ,  AND  SA2  OE  STIFFNESS  * 
C  MATRIX  , 

CALL  AMTRX4 
WRITE(6,501)AG  ,  , 

501  FORMAT?//  '  AG:  ’.(12F4.1)) 

WRITE?© ,S04)BG  , 

504  FORMAT?//  BG:  ’,(12F4.I)) 

WRITE?© , l002 )GA1 

1002  FORMAT?//  GAl  ,  /  ,  ( 3X ,  6F7 . 3 )  ) 

WRITE? © . l012 )SA1 

1012  FORMAT?//  ,  SAlt ,  /  ,  (3X , 6F7 . 3 )  ) 

,  WRITE?6.l004)GBl 
1004  FORMAT?  /  7/  GB1  ’  ,  /  ,  (  3X ,  6F7 . 3 ) ) 

,  WRITE X © , 10 14 ) SB 1 
low. .  FORMAfll/ •:  'SBI '  /  ,  (3X, 6F7 . 3 ) ) 


1002 

1012 

1004 


1014 
C  L 


FORMAl 

WRITE? 

FORMAl 

WRITE? 

FORMAl 

WRITE? 

FORMAl 

WRITE? 

FORMA! 

WRITE? 

FORMAl 


’ , (12F4. 1) ) 
' , ( 12F4 . 1 ) ) 


10 14) SB 1 
/,  f  SBl’y 

Sector  wi 


LOAD  BORDER  VECTOR  W 

CA^  BORDER fWl,V) 
WRIT??6.510)V 

li=4*nl6ng 


WRITE /6?520?fWl(L) ,L=1.L1) 

ffifiS? \ti2, V 

FORMAT ?3X  5F8. 2) 

FORMAT?///  ,V:  4 ,/.  (4X.6F10.5)) 
CALL  MULtUWl.y^Al.SBi)  4 
WRITEC6 ,520) ?^1?L) ,t=l,Ll) 
WRITE?6  521)?W1?L)  L=L2,L3) 
WRITE?6 ISlOjv  / 

CALL  BOADERfwi.V ) 4 
WRITE?6,520)(Wi?L) ,L=1.L1) 

WRITE ?6j 52 lj?Wl?L) ’L=l2,l3) 

S^^t5t4W'v,s*1'GBt) 

WRITEtl  )§20KW1(L) ,L-1,L1) 
WRITEi$.s5o)P 

FORMAT?//  P:  '  ,  /  ,  (6X,  3F12 .4)  ) 


/;?3X,6F8:2)) 

/.(4X.6F10.5)) 


531  FOR] 
C 

C  FORM  U 
C 


READ(5  *)D 
WRITeTS ,531  )D 

format!/  ,  D: 


.3F12.4) 


VINT*P 


LONGM=NLONG- 1 
LONGMM=  NLONG - 2 
NLATM=NLAT- 1 
DO  30  L= 1 , NLATM 

jpi=!l-iHnlatm 

KUl=(L-l1*(NLONG-2)-l 
DO  29  K=2,LONGM 
TEMP=0 . 

DO  2§  J=l, NLATM 

jv=j*nlon6+k 

JP=JP1+J 

28  TEMP=TEMP+V (JV)*P ( JP ) 

29  U  7KU 1 + K )  = TEMP 

30  CONTINUE 
WRITE(6,532)U 

532  FORMATl?  ,Y/U:  '  /  .  ( 6X.4F12 . 4) ) 

CALL  FFf6b(U , D,GAi, 9a1) 

WRITE (6, 532)6 

C  PUT  FINAL  RESULT  IN  V 
C 

DO  40  L=l, NLATM 
DO  39  K=l,LONGMM 
TEMP=0 . 

DO  38  J  =  1 , NLATM 

3  8  TEM? = TEMP  +  P(CJ- 1)*NLATM+L )  *U  (  (  J  - 1 ) *LONGMM+K ) 

39  V(L*NLONG+kU)=TEMP 

40  CONTINUE 
WRITE (6, 5 10 )V 
STOP 

END 

C 

c*********************************************************** 

************************************ 

c 

C  THIS  SUBROUTINE  CLEARS  THE  BORDER  VECTOR  W1  AND 
C  SUBSTITUTES  THE  BOUNDARY  VALUES  FROM  V. 

IMPLICIT  REAL*8 (A-H.O-Z) 

COMMON/ CM1A/NLAT, NLONG 
DIMENSION  Wl(ZQ)  V(ZP) 

LATX=NLAT+ 1 
NC=4*NL0NG 
NR=4*LATX 
NB=NC+NR 

,  DO  4  J=1,NB 

4  W1(J)=0. 

DO  6  J  = 1 .NLONG 
6  Wl(J)=vm 

DO  8  J=l, NLONG 
L=3*NLQNd+ J 
K=NLAT*NLONG+J 
8  W1(L}=V(K) 

DO  10  J=1,LATX 
L=NC+J 

K= ( J- 1 )*NLONG+ 

W1(L)=V(K) 

DO  12  J  *  1  j. LATX 
L=NC+3*LAtx+J 
K=J*NLONG 


RETUi 


Q  **********************TWf  ********************************* 

c  ****§5S52liJJ§5**^51^&************************************ 

C  THIS  SUBROUTINE  FORMS  THE  MATRICES  GAl .  GB1 .  SA1,  AND  SB1 
C  THAT  ARE  FACTORS  IN  THE  TENSOR  PRODUCT $  USE6  TO  #ORM  THE 
C  COEFFICIENT  MATRIX  ("STIFFNESS"  MATRIX)  FOR  THE  POISSON 


c 

c 

c 


c 

c 

c 


EQUATION.  ALL  OF  THESE  MATRICES  ARE  SYMMETRIC  AND 
TRIDIAGONAL. 


1GA1T3"'NL0NG- 


LATX=NLAT+1 
LONGM=NLONG- 1 

FIND  BG  =  (ELEMENT  HEIGHT)/ 6. 
NM=NLONG- 1 


SAl(ZK) ,GBI 

gbi(2*n1at- 


£ZL) , SB1 ( ZL) 


SUBROUTINE  MULT1(W1,V,A,B) 

SUBROUTINE  PREMULTIPLIES  W1  MATRIX  BY  TRUNCATED  A  MATRIX 
(FIRST  AND  LAST  ROWS  OMITTED),  POSTMULTIPLIES  PRODUCT  BY 
TRUNCATED  B  MATRIX  (FIRST  AND  LAST  COLUMNS  OMITTED),  AND 
SUBTRACTS  INTERIOR  ELEMENTS  OF  W1  FROM  CORRESPONDING 

Q**§5c53S^5§'*95*X'i************'*  ■*■***•*■**•*■*****•****■*••*  •****■**■***** 


c 

c 

c 

c 

c 

c 

c 

c 


"oiA-n.u- 

lat,nl6ng 

zqj;vTzpj,a(1),b(1) 


IMPLICIT  REAL*8 (A-H.O-Z) 
COMMON/ CM1A/NI  *”  A”~ 
DIMENSION  Wl(2 
LATX=NLAT+ 1 
LONGM=NLQNG-l 
FCU=W1(1) 
RCU=Wl(3^NLONG+l) 
L1=4*NL0NG 

L3=LlU*LATX 
DO  2  J=2 .LONGM 
K=2*(J-l} 
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L=3*HL0NG+J 

FCC=a7k+1)*FCU+A(K)*W1( J)+A(K+3)*W1( J+l) 
RCC=A(K+i)*RCU+A(K)*Wl(L)+A(K+3)*Wl(L+l) 
FCU=WI(J) 

RCU=V1IL) 

wiTj)=fcc 

W1[L)=RCC 
LAS  T  A = 2  N LONG  - 1 
DO  4  J=2,NLAT 
L=4*NL0N6+J 


ll=l+latx 

K=L+3*LATX 
KK=K-LATX  , 
W1(LLT=A(3)*W1 
W1(KK)=A(LASTA 
WRITE ( 6 ,  $20 ) ?W 

r  tn  rmn  >  /•  '  r  A  1  \  >  »  i 


WRITER 

FORMATI 

FORMAT 

LASTB=i 


=4*NL0NG 


521)(Wia 
>//  INTE1 
JX  6F8.2) 
^lAtx- i 


!l=1.L1) 

,L=L2,LJ ) 

MEDIATE  RESULT,  Wl* , / , (3X, 5F8 . 2) 


gl CNLONGV 2£=B f 3 ) *W1 ( 2 ) +B ( 2 ) *W1 (NC+LATX+ 2 ) +B ( 5 ) 

Wl (2*NLONG- 1} ?B ( 3) *W1 (NLONG- 1 ) +B ( 2 ) *W1 (NC+ 2*LATX+ 2 ) 
+B  5 )*W1 (NC+2*LATX+  3  J  . 

Wl  2*NLONG+2)=B?LASTB-2)*yi{NC+2*LATX-2)+B(LASTB-3) 
*W;,TnC^2*LATX- 1 )  fB(LASTB  )*W1  f  3*NL0NG  +  2) 

Wl  3*NL0NG- 11=B (LASTB- 2pfVlfNC+  3“LATX- 2 ) +B (LASTB-3 ) 


;nc+2*latx- 1 ) +b(last] 

Wl(3*NLONG-l>B(LASTB-2)*Wl(NC+3*LAT] 
L*W1CNC+  3*LATX- 1 ) +  B (LASTB ) *WI ( NC- 1 ) 
J2=NLON G-2 

Wlf  J + S^NLolcJ = B  LAS  £  B  )  *W1  (  J  ♦  3  *NL0NG  ) 


!*yi(NC+2*LAT 
h*Wi(3*NL0NG 
P-W1{N0  3*LAT 


URL=Wl(NC+LAtX+2) 

BRL=Wl(NC+2*LATX+2) 

J2=LATX-2 
DO  8  J=3 , J2 
ND=2*(J-i) 

NCU=NC+LATX+J 

URC=B(ND+ 1 ) *URL+B(ND } *W1 (NCU ) +B(ND+ 3 )*W1 (NCU+ ! 
BRC=B (ND+ l1*BRL+B (ND ) *W1 (NCU+LATX ) +B (ND+  3 ) 
l*Wl(NCy+LAfX+l7  v 


URL=W1(NCU) 
BRL=W1( NCU+LATX) 

wi7ncu)=urc 

Wl (NCU+LATX) =BRC 
W1(NC+LATX+2)=W1 
Wl(NC+2*LATX-l)= 
WllNC+2*LATX+2 )= 
W1(NC+3*LATX-1)= 

CORRECT  V 

NLATM=NLAT- 1 
DO  10  J = 3 , NLATM 
L=NC+LATX+ J 
K=(J-l|*NL0NG+2 

L=LATX+L 
K=J*NLQNG-1  .  . 


ONG+2) 
2*NLONG+2 
2*NL0NG- 1 
3,VNL0NG-  1 


K=J*NL0NG-1  #  x 

DO  12  J=2,J2 

V(S=V?lJ-W1(L) 

L=2*NLONG+y 

K=?NLAT-1)*NL0NG+J 

VXk)=V(K)-W1(L) 

RETURN 


C***********************VWf********************************* 

SUBROUTINE  FFFDB (X , E , GA , SA ) 

C  THIS  SUBROUTINE  SOLVES  A  SUCCESSION  OF  ONE - DIMEN  SION SAL 
C  PROBLEMS.  THE  RELEVANT  COEFFICIENT  MATRIX  C  IS  FIRST 
C  FORMED,  THEN  FACTORED,  FOLLOWED  BY  FORWARD  REDUCTION, 
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